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The Airy distribution function describes the probability distribution of the area under a Brownian 
excursion over a unit interval. Surprisingly, this function has appeared in a number of seemingly 
unrelated problems, mostly in computer science and graph theory. In this paper, we show that this 
distribution function also appears in a rather well studied physical system, namely the fluctuating 
interfaces. We present an exact solution for the distribution P(/i m ,L) of the maximal height h m 
(measured with respect to the average spatial height) in the steady state of a fluctuating interface in 
a one dimensional system of size L with both periodic and free boundary conditions. For the periodic 
case, we show that P(h m , L) = L" 1 ^ 2 f (/i m L -1 / 2 ) for all L > where the function f(x) is the Airy 
distribution function. This result is valid for both the Edwards- Wilkinson and the Kardar-Parisi- 
Zhang interfaces. For the free boundary case, the same scaling holds P(h m , L) = L~ 1/2 F (h m L~ 1/2 ) , 
but the scaling function F(x) is different from that of the periodic case. We compute this scaling 
function explicitly for the Edwards- Wilkinson interface and call it the J 1 - Airy distribution function. 
Numerical simulations are in excellent agreement with our analytical results. Our results provide 
a rather rare exactly solvable case for the distribution of extremum of a set of strongly correlated 
random variables. Some of these results were announced in a recent Letter [ S.N. Majumdar and A. 
Comtet, Phys. Rev. Lett., 92, 225501 (2004) ]. 

PACS numbers: 02.50.-r, 89.75.Hc, 89.20.Ff 

I. INTRODUCTION 

A Brownian excursion x(r) is a conditioned one dimensional Brownian motion over the time interval < r < T 
such that its path starts and ends at the origin x(0) = x(T) = 0, but is constrained to stay positive in between 
(see Fig. 2). The area under the excursion, A = f x(t)(1t, is clearly a random variable taking a different value 
for each realization of the excursion. A natural question that the mathematicians have studied quite extensively 
over the past two decades is: what is the probability distribution P(A, T) of the area under a Brownian excursion 
over the interval [0, T]? Since the typical lateral displacement of the excursion at time r scales as yfr, it follows 
that one can trivially rescale x(t) = ^/ti/(t/T) where y(u) represents a Brownian excursion over the unit interval 
u e [0, 1]. Hence, the area A = J x(r)dT ee T 3 / 2 Jg 1 y(u)du scales as T 3 / 2 . Thus the area distribution has a scaling 
form, P(A,T) = T- 3 / 2 f (A/T 3 / 2 ). The normalization condition / °° P(A, T)dA = 1 demands a prefactor T~ 3 / 2 and 
also the conditions: f(x) > for all x and / °° f(x)dx = 1. One then interprets the scaling function f(x) as the 
distribution of the area under the Brownian excursion y(u) over a unit interval u £ [0, 1]. The function f{x), or rather 
its Laplace transform, was first computed analytically by Darling [1] and Louchard [2], 

f(s) = / f{x)e- sx dx = sV2^ V e -^ s2/32 ~ 1/3 , (1) 
Jo k=i 

where a^'s are the magnitudes of the zeros of the standard Airy function Ai{z) on the negative real axis. For example, 
u\ = 2.3381 ...,a2 = 4.0879 . . ., = 5.5205 . . . etc. Since the expression of f(x) involves the zeros of Airy function, 
the function f(x) has been named the Airy distribution function, which should not be confused with the Airy function 
Ai(x) itself. Even though the Eq. (1) provides a formally exact expression of the Laplace transform, it turns out the 
calculation of the moments M n = J °° x n f(x)dx is highly nontrivial and they can be determined only recursively [3] 
(see Section II). Takacs was also able to formally invert the Laplace transform in Eq. (1) to obtain [3], 

A*) = ^jte- bk/x2 bl /3 U(~5/6A/3,b k /x 2 ), (2) 
fc=i 

where = 2a|/27 and U(a, b, z) is the confluent hypergeometric function [4]. The function f(x) has the asymptotic 
tails [3,5], 
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f(x) ~ x- 5 e - 2a *' 27x2 as x^O 

f(x) ~ e~ 6x2 as x -> oo. (3) 
A plot of this function, obtained by evaluating the sum in Eq. (2) using the Mathematica, is provided in Fig. 1. 
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FIG. 1. A Mathematica plot of the Airy distribution function f(x) in Eq. (2). 

So, why would anyone care about such a complicated function? The reason behind the sustained interest and study 
[3,5-8] of this function f(x) seems to be the fact that it keeps resurfacing in a number of seemingly unrelated problems, 
mostly in computer science and graph theory [for a list of such problems see [6]]. For example, the function f(x) 
describes the distribution of the cost of the construction of a linear table for data storage using the linear probing with 
random hashing algorithm [6]. The function f(x) also describes the distribution of the total path length in Catalan 
trees [3]. The generating function for the number of inversions in trees involves the Airy distribution function f(x) [9]. 
Also, the moments M n 's of the function f(x) appear in the enumeration of the connected components in a random 
graph [10,11]. Recently, it has been conjectured and subsequently tested numerically that the asymptotic distribution 
of the area of two dimensional self-avoiding polygons is also given by the Airy distribution function f(x) [12]. Besides, 
numerical evidence suggests that the area enclosed by the outer boundary of planar random loops is also distributed 
according to the Airy distribution function f(x) [13]. Given the widespread occurrence of this function fix) in various 
combinatorial problems, it is natural to ask whether this function also appears in a more realistic physical system. In 
a recent Letter [14], we have demonstrated that indeed this Airy distribution function also appears in a rather well 
studied physical system, namely the one dimensional fluctuating interfaces. 

Fluctuating interfaces have been widely studied over the last two decades as they appear in a variety of physical 
systems such as growing crystals, molecular beam epitaxy, fluctuating steps on metals and growing bacterial colonics 
[15,16]. The past theoretical efforts mostly concentrated in characterizing the universal scaling properties of the 
surface roughness measured via the average width of the surface. This average roughness, however, fails to measure 
the possible extreme fluctuations which, though rare, may play an important role under many situations. For example, 
extreme fluctuations play a significant role in batteries where a short circuit may occur if the highest point of the 
metal surface on one electrode reaches the one on the opposite electrode. There is also a compelling theoretical 
motivation for studying the extreme fluctuations. While the extreme value statistics is well understood for a set of 
independent random variables [17], it becomes nontrivial when the random variables are correlated. The extreme value 
problems associated with correlated variables have appeared recently in a variety of problems ranging from disordered 
systems [18] to a number of computer science problems that include the growing search trees [19] and the networks 
[20] . A fluctuating interface is an ideal candidate for studying the extreme statistics of correlated variables since the 
heights at two different points on the interface are correlated. A knowledge of the distribution of their maximum (or 
minimum) may possibly provide important insights into this important general class of extreme value problems of 
correlated variables. 

Motivated by these observations, Raychaudhury et. al. [21] recently proposed to study the statistics of the global 
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maximal relative height (MRH) (measured with respect to the spatially averaged growing height) of a fluctuating 
(f + f )-dimensional interface. By (1 + l)-dimcnsional interface one means an interface characterized by a scalar height 
function H(x, t) growing over a linear substrate. For a substrate of finite size L, the system reaches a stationary 
state in the long time limit. In this stationary state, the MRH h m is expected to scale as the surface roughness, 
h m ~ L a for large L, where a is the roughness exponent. This suggests, quite generically, a scaling form for the 
normalized probability distribution of the MRH, P(h m ,L) <~ L~ a f\ (h m /L a ) where fi{x) is a scaling function. This 
was demonstrated numerically in Ref. [21] for a one dimensional lattice model belonging to the Edwards- Wilkinson 
(EW) universality class [22]. Further, it was argued that the scaling function fi(x) is sensitive to the boundary 
conditions [21]. 

In our previous Letter [14], using path integral techniques, we presented an exact solution of the stationary MRH 
distribution P(h m ,L) for the (1 + l)-dimensional EW model, with both periodic and free boundary conditions. For 
the periodic case, we showed that P(h m ,L) = L^ 1 / 2 f {h m L~ x l 2 } for all L > where the scaling function f(x) is 
precisely the Airy distribution function defined in Eq. (1). For the free boundary case, the same scaling was found 
to hold though the scaling function was different from that of the periodic case. The purpose of the present paper is 
twofold: (i) to provide detailed derivations of various results that were just announced in Ref. [14] and (ii) to extend 
our results to other systems such as the Kardar-Parisi-Zhang (KPZ) interface [23]. Our main results, along with a 
layout of the paper, are summarized below. 

1. In Section II- A, we provide a new physical derivation, using methods of path integral, of the distribution of the 
area under a Brownian excursion. This result is used later in Sec. IV-A in the context of the MRH distribution of 
an EW interface with periodic boundary condition. In Section II-B, we derive the distribution of the area under 
a Brownian meander (a Brownian meander is a one dimension Brownian motion x(r) over the time interval 
< r < T that is pinned to the origin at one end x(0) = 0, but is free at the other end at r = T, and is 
constrained to stay positive in between). These results for the Brownian meander will be used later to discuss 
the MRH distribution of an EW interface with free boundary condition in Sec. IV-B. The Section II can be 
read independently of the rest of the paper and the results here are interesting by themselves. 

2. Section III includes a detailed derivation of various correlation functions in the EW model in two opposite 
regimes: (a) growing regime where the correlation length of the interface £(t) <~ t 1 ^ (z being the dynamical 
exponent) is less than the system size, £(t) << L and (b) stationary regime when >> L and thus the 
joint probablity distribution of height fluctuations become time independent. An explicit form of this joint 
distribution in the stationary regime is derived for both periodic and free boundary conditions. 

3. In Section IV, we derive analytically the stationary MRH distribution in (1 + l)-dimensional EW model. Section 
IV-A discusses the periodic boundary condition case, where the scaled MRH distribution is shown to be exactly 
the Airy distribution function defined in Eq. (1). The moments of this distribution are also computed exactly. 
The analytical results are then compared to the numerical results obtained from the direct numerical integration 
of the EW equation and an excellent agreement is found. In Section IV-B, we derive the MRH distribution for 
the EW interface with free boundary condition and show that it is related to the area under two independent 
Brownian meanders. We call the associated scaling function F-Airy distribution function where F refers to 
the free boundary condition. We also compute the moments and the asymptotic properties of this distribution 
analytically. As in the periodic case, we find excellent agreement between the analytical and the numerical 
results. 

4. In Section V, we discuss the (l + l)-dimcnsional KPZ equation and show that for the periodic boundary condition, 
the stationary MRH distribution is again given by the Airy distribution function in Eq. (1). This is further 
confirmed by direct numerical integration of the (1 + l)-dimensional KPZ equation. 

5. Section VI discusses the MRH distribution in the growing regime t « L z . Using the well known theory 
of the extreme value statistics of independent random variables, we show that the appropriately scaled MRH 
distribution is universal and is given by the Gumbel function. 

6. We conclude in Section VII with a summary, experimental consequences of our results and some open questions. 

II. THE AREA UNDER A BROWNIAN EXCURSION AND A BROWNIAN MEANDER 

In this section, we provide a new derivation, using path integral methods, of the distribution of the area under a 
Brownian excursion and a Brownian meander. These results were originally derived by mathematicians [1-3,24,25] 
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using various probabilistic methods. The path integral method presented here, though perhaps not rigorous in the 
strict mathematical sense, does provide a simpler physical derivation. Thus the results of this section are interesting 
by themselves and can be read independently of the rest of the paper that discusses the fluctuating interfaces. This 
section will also serve as an introduction to the path integral method to be used in later sections to study the maximal 
height fluctuations in interfaces. 



A. Brownian Excursion 



A Brownian excursion is a one dimensional Brownian motion over the time interval < r < T that starts and 
ends at the origin x(0) = x(T) = 0, but is constrained to stay positive in between. We are interested in computing 
the probability distribution P(A, T) of the area A = J Q x(r)dr under the excursion. The Brownian excursion, as 
described above, is however a bit ill defined for a Brownian walk in continuous space and time. For such a walk, it 
is well known that if the walker crosses zero once, it recrosses the zero infinitely many times immediately after the 
first crossing. Thus, for a continuous space-time Brownian motion, it is impossible to enforce the constraint x(0) = 
and simultaneously forcing it to stay positive immediately after. There are two ways to go around this problem. The 
cleanest way to study an excursion is to consider a discrete time random walk moving on a discrete one dimensional 
lattice and then sample all configurations that start at the origin and come back to the origin for the first time after 
an even T = 2n number of steps. Such paths are called Dyck paths [3]. One then defines the area under such a path 
as A — J2i2i x i where Xi is the position of the walker at step i. By appropriately taking the continuum limit, one 
would then arrive at the Brownian excursion [26]. This method, though conceptually clear, is somewhat cumbersome 
mathematically. 

Here we devise a second method that turns out to be more amenable to path integral techniques. We consider the 
Brownian motion in continuous space and time, but introduce a small cut-off e by hand as shown in Fig. 2. More 
precisely, our Brownian motion starts at x(0) = e and comes back after time T to x{T) — e, without crossing the 
origin in between. We will first derive the probability density P(A, T, e) of the area A — J Q x{t)<1t for a fixed e 
and then finally take the limit e — > 0. We will show that the limiting distribution exists and is precisely the Airy 
distribution function defined in Eq. (1). 




FIG. 2. A Brownian excursion over the time interval < r < T starting at x(0) 
positive in between. 



e and ending at x(T) = e and staying 



The probability measure associated with an unconstrained Brownian path x(t) over the time interval < t < T is 
proportional to exp — i (^f ) 2 dr . In addition, we need to incorporate the constraint that it stays positive between 

and T which can be achieved by multiplying the above measure with the indicator function IIt=o ^ l x ( T )] which is 

1 if the path stays positive in < r < T and zero otherwise. The distribution P(A,T) of the area A = J Q T x(r)dT 
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under this Brownian excursion can then be expressed as a path integral 

rx(T)=e rT T ( rT 

Vx(T)e--J« dT ^' d ^ TT [ x (t)] S / 

/x(0)=e T=0 \Jo 



P{A,T) = ^- Vx{T)e—^o dT{dx,dT) X\6[x{T)]8\ x(r)dr - A . (4) 

Jo;(0)=e „_ n Wo / 



Note that we have suppressed the e dependence of P(A, T) for brevity. The normalization of the distribution, 
J Q °° P(A, T)dA — 1, is ensured by the following definition of the partition function Ze of the Brownian excursion 



r x(T)=e f T T 

Z E = Vx{T)e-*J a dT{dx/dT) \[e[x(T)]. (5) 

Jx{0)=e r=0 



All paths inside the path integrals in Eqs. (4) and (5) propagate from their initial value x(0) = e at r = to their 
final value x(0) — e at r = T. 

We first evaluate the partition function .Ze. This is easy since one can identify the quantity inside the exponential 
in Eq. (5) as the action corresponding to a single particle quantum Hamiltonian, Hq = + Vo(x), where the 

potential Vo(x) = for x > and Vo(x) — oo for x < 0. The infinite potential for x < ensures that the path 
never crosses zero and thus takes care of the indicator function I1t=o ^ i x ( T )] m ^q. (5). Using the standard bra-ket 
notation, the partition function Ze is then simply the propagator 

Z E =< e\e- HaT \e > . (6) 

It is easy to see that the Hamiltonian H has continuous eigenvalues Eq = k 2 /2 labelled by k > and the corresponding 
eigenfunctions that vanish at the origin arc given by, ipk(x) = ^J^sin(kx). Expanding the propagator in the energy 
eigenbasis, one gets 

POO -I 

Ze = 1 dk\Me)\ 2 e-* T ' 2 = ^=(l-e-"' T ). (7) 

Note that Ze is just the probability that a Brownain motion goes from x(0) = e to x(T) = e in time T without 
crossing the origin, and hence can also be computed by the standard image method which yields the same answer [27] 
as in Eq. (7). 

We now turn to the cvalution of P(A, T) in Eq. (4). Taking the Laplace transform, P{\, T) = / °° P(A, T)e- XA dA 
where A > 0, in Eq. (4) gives, 

i />x(T)—e r T ri 2 l — 

P( A , T) = J- / Vx(t) e~ Jo dT TT e [x(r)] . (8) 

J 0)=* r=0 

The numerator in Eq. (8) can then be identified as the propagator < e\e~ HlT \e >, with the Hamiltonian Hi = 
~\~Ix I + Vi( x ) wnere V\(x) is a triangular potential: V\ (x) = Ax for x > and V\ (x) = oo for x < 0. The later 
condition again takes care of the fact that the paths do not cross zero. The Hamiltonian Hx has only bound states 
and hence discrete eigenvalues. Solving the corresponding Schrodinger equation, one finds that the wavefunction (up 
to a normalization constant), that vanishes as x — > oo, is simply given by, iPe(x) = Ai [(2A) 1 / 3 (x - E/Xj] , where 
Ai(z) is the standard Airy function [4]. The condition that the wavefunction should vanish at x = determines the 
discrete eigenvalues, E& = afeA 2 / 3 2 -1 / 3 for k = 1,2, . . ., where a^s are the magnitude of the zeros of Ai(z) on the 
negative real axis. For example, one has ot\ — 2.3381 . . ., a 2 = 4.0879 . . ., as = 5.5205 . . . etc. Thus the normalized 
eigenfunction is given by 

M[(2X)^x-a k ] 

i>k{x) = ; J ^=- (9) 



y/f™Ai? [{2\)V*y-a k ]dv 



The integral J °° Ai 2 [(2A) 1 / 3 y — a^] dy = (2A) 1/13 J X ' ak Ai 2 (z)dz can be further simplified by using the identity, 

J x ak Ai 2 (z)dz = [Ai'(-a k )} 2 [28] where Ai'(z) = dAi(z)/dz. Expanding the propagator < t\e" HlT \e > into the 
energy eigenbasis we get 
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< e\e- H ^ T \e >= £ \M^\ 2 e~ X 2 akT - (10) 
fe=i 



Using the exact eigenfunctions from Eq. (9) and then substituting all the expressions in Eq. (8) we get 

We are now ready to take the e — > limit. Expanding both the numerator and the denominator in Eq. (11) in a 
Taylor series in e and taking the e — > limit, we get a simplified result, 

oo 

P(A,T) = V2^(AT 3 / 2 )^ e - 2 " 1/3 ^( AT3/2 ) 2/3 . (12) 
fc=i 

Note that the right hand side of Eq. (12) is a function of only one scaling combination s = AT 3 / 2 . This demands that 
the distribution P(A, T) must have the scaling form P(A, T) = T~ 3 / 2 / (A T~ 3 / 2 ) , so that its Laplace transform is a 

function of only the combination s — AT 3 / 2 , P(X, T) = J °° f{x)e~ xr3 x dx. Comparing this to the right hand side of 
Eq. (12) we arrive at the final result 



r oo 

/(«)=/ f{x)e- sx dx = sV2^y^e- a > 
J ° fc=i 



s 2/3 2 -l/3 



(13) 



Moments of the area distribution: Using the scaling form, P(A, T) = T 3 / 2 / (A T 3 / 2 ) , one finds that the 
moments A n = J °° A n P(A,T)dA = M n T Zn / 2 , where M n — / °° f(x) x n dx arc the moments of the area under the 
excursion over a unit interval. The extraction of the moments M n explicitly from the Laplace transform in Eq. (13) 
is highly nontrivial. However, starting from Eq. (13), Takacs found a recursive method to compute the moments [3]. 
We summarize the main results here without the details. One first defines a set of new variables K n via the relation 

M n ^^-^^^K n , (14) 

where T{x) is the Gamma function. The variables K n : s subsequently satisfy a nonlinear recurrence 

3 — 4 ™ _1 

K n - -^—K n ^ + KjKn-j, (15) 

starting with K a = —1/2. This produces recursively the exact values of the moments M n for any n. For example, the 
first few values are, 

*, 1 /tt ,.5 15 fn %r 221 

M = 1, M 1 = -J-, M 2 = —, M 3 = — J-, M 4 = , ... (16) 

U 2 V 2 ' 12' 3 64 V 2' 1008' V ' 

These moments will be used later in Sec. IV-A in the context of the maximal height fluctuations in interfaces with 
periodic boundary conditions. 



B. Brownian Meander 



A Brownian meander is a path of a one dimensional Brownian motion over the time interval < t < T that starts 
at the origin x(0) = and stays positive upto r = T. The difference between an excursion and a meander is that in 
the former case, the path comes back to the origin at the end of the interval, x(T) = 0. In the case of the meander, 
the Brownian walker is free to arrive at any final position as long as the final position is positive. As in the case of 
the Brownian excursion in the previous subsection, we consider the Brownian meander in continuous space and time, 
but introduce a small cut-off in the initial position x(0) = e, and eventually take the e — > limit. In Fig. 3, we show 
a Brownian meander that finally arrives at the position b > at time r = T, starting at x(0) = e. Contrary to the 
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Brownian excursion where the final position is same as the initial position i.e., x(T) = b = e, for the meander we need 
to integrate the paths over all possible final positions b > 0. 




FIG. 3. A Brownian meander over the time interval < r < T starting at a;(0) = e and ending at x(T) — b and staying 
positive in between. 

As before, we are interested in computing the probability distribution P(A, T) of the area A — J„ T x{r)dr under a 
meander. We suppress the e and b dependence of P(A, T) for brevity. Following the same arguments as in Sec. H-A, 
one can easily write down an expression for the Laplace transform, P(A, T) = J °° P{A, T)e~ XA dA, as a ratio of two 
propagators, 



P(X,T) 



JT db < b \ e 



—H\T 



e > 



f™db<b\e-"° T \e> ' 



(17) 



where the Hamiltonians H and Hi are the same as in the previous subsection. 

The denominator, Z M = J °° db < b\e~ H ° T \e >, is the partition function of the Brownian meander and can be easily 
evaluated by decomposing the propagator in the energy eigenbasis of H . We get 

Z M = erf(e/\/2T), (18) 

where erf (a;) = f£ e~ u ' ' du is the standard Error function. Note that Zm is simply the probability that a Brownian 
walker does not return to the origin upto time T starting at the initial position e and hence, can be computed also 
via the method of images [27]. 

The propagator in the numerator in Eq. (17) can be evaluated by expanding it in the energy eigenbasis of H\, 
which were detailed in sec. II-A. We do not repeat the calculations here and just mention the final result. We get 



P(X,T) = 



E 



Ai[{2\)^e-a k ] J~ Ai(z)dz 



_ A 2/3 2 -l/3 T 



(19) 



crf(e/V2T) ^ M'\-a k ) 
where a^s are the magnitude of the zeros of the Airy function. Taking the e — > limit we obtain a simpler expression 



3/2N2/3 



P(A,T) = V?F2- 1 /6 (AT 3/ 2) i/3 Y^B{a k )e- 2 ' a ^ XT ^ 



fc=i 



where 



B(a k ) = 



Ai'(-a k ) 



(20) 



(21) 
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One should compare this result with its analogue in the case of the Brownian excursion in Eq. (12). As in the case of 
the excursion, the right hand side of Eq. (20) is a function of only one scaling combination, s = AT 3 / 2 . This indicates 
that P{A,T) = T~ 3 / 2 g (A T -3 / 2 ) , where the Laplace transform of the scaling function g(x) is obtained from Eq. 
(20), 

g(s)= I g(x)e- sx dx = V^^ 1/e s 1/3 y B(a k )e- a ^" 2 ~ L '\ (22) 



p OO _ _ 

= g(x)e- sx dx = V^2~ 1/6 s 1/3 ^2B(a k )e 
Jo fc =i 



To make contact with the existing results in the probability literature, we first note that the Laplace transform, 
P(X,T) = P(A,T)e~ XA dA can be regarded as the expectation value E[e~ XA ] with respect to the probability 
distribution P(A, T). Next, in terms of a new scaling variable t = 2- 1 / 3 X 2 / 3 T, Eq. (20) can be rewritten in a simpler 
form, 



E 



ntJ2 B (<*k)e- akt , (23) 



fe=i 



where we have defined the scaled area a = A/T 3 / 2 . Note that, as in the case of an excursion, the interval length T 
only appears as a trivial scale factor and hence one can interpret a as the area under a Brownian meander over the 
unit interval [0, 1]. Dividing both sides of Eq. (23) by \/nt and taking a further Laplace transform with respect to t, 
we get 



-y/2t 



\/wt L J (a* + u) Ai(u) 

The second equality in Eq. (24) can be proved by replacing the sum by a contour integration in the complex plane, 
a derivation of which is given in Appendix-A. This result in Eq. (24) is identical to an expression derived by Perman 
and Wcllncr [24] using probabilistic methods. Here we have provided a simpler physical derivation using path integral 
techniques. 

Moments of the area distribution: The scaling form P(A,T) = T~ 3 I 2 g (AT -3 / 2 ) indicates that the moments 
of the area under a Brownian meander are given as A n = J °° A n P{A, T)dA = a n T 3n l 2 , where a n = J^g{x)x n dx 
are the moments of the area under a Brownian meander over the unit interval. As in the case of Brownian excursion, 
the extraction of the moments a n 's explicitly from the Laplace transform in Eq. (22) is not easy. However, starting 
from the alternative expression of the Laplace transform in Eq. (24), Perman and Wellner were able to obtain the 
moments a„'s recursively [24]. Here we briefly summarize their results without the details. We first define a set of 
new variables i?„'s via the relation 

«« = V^2-"/ 2 ^^i?„. (25) 
Next, the RnS were found to satisfy the following recursion relations [24] for all n > 1, 

n 
3 = 1 

3 

A» = 7n + ^(2n-l))8 n _i 

(36)^" r(3n+l/2) 
ln ~ r(n+l) r(n + l/2) " { ' 

Using Eqs. (26) and (25), the moments a„'s can be subsequently calculated recursively. The first few values are 

3 pi: 59 465 pi: 

00 = 1, "1 = 4^2' a2 =60' a3= 512V2' - (27) 

We will use these results later in the context of the calculation of the moments of the maximal height fluctuations in 
interfaces with free boundary conditions. 
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III. FLUCTUATING EDWARDS- WILKINSON INTERFACE 



In this section, we will discuss some general properties of the dynamics and the stationary state of a fluctuating (1 + 
l)-dimensional interface. We consider here the simplest model of a fluctuating one dimensional interface characterised 
by a single valued height function H (x, t) which is growing on a linear substrate of size L according to the celebrated 
Edwards- Wilkinson equation [22] 

^ = (28) 

where r](x,t) is a Gaussian white noise with zero mean and a correlator, (r](x,t)r](x't')) = 25(x — x')8(t — t'). The 
equation (28) has a soft (zero wave vector) mode since the spatially averaged height H(x, t) — J Q L H (x, t)dx/L diffuses 
with time (typical value growing as \/t/L) even in a finite system of size L. Thus, the joint height distribution P[{H}, t] 
of the system will never reach a time independent stationary state. Hence, it is useful to subtract this zero mode from 
the height and define the relative height, h(x, t) — H(x, t) — H{x, t), whose joint distribution P[{h}, t] will eventually 
reach a stationary state in the long time limit in a finite system. Note that the average height H (x, t) is over a given 
sample, and not the statistical average over the noise r\ in Eq. (28). So, from now on, we will always consider the 
relative height h(x, t) as our basic variables. Note that, by definition, 

L 

h(x, t)dx = 0. (29) 

We will see later that this constraint of zero total area under the relative height h plays an important role in determining 
the maximal relative height (MRH) distribution. All the other nonzero modes of h evolve identically as those of the 
actual height H. 

We will see shortly that there are two regimes of the evolution of the interface in a finite system. Let us assume 
that we start from a flat initial configuration where the heights at different points on the substrate are completely 
uncorrelated from each other. As time t grows, the heights at different points get correlated and the system is 
characterized by a single growing correlation length £(t) ~ i 1 / 2 , where z is the dynamical exponent (z = 2 for the 
EW interface). As long as this growing correlation length is much smaller than the system size, £(t) ~ << L, the 
interface is in the growing regime. In this regime, the system does not feel the boundary and hence is insensitive to the 
boundary conditions. When the correlation length becomes comparable or exceeds the system size, 
the system crosses over to a stationary regime where the joint distribution of relative heights P[{h},t] becomes time 
independent. In this regime, the system is strongly sensitive to the boundary conditions. The crossover time between 
the two regimes scales as t c <~ L z . 

The two subsections below deal with two different boundary conditions, respectively the periodic and the free 
boundary condition. Our goal in these subsections would be to compute the height-height correlation function as well 
as the stationary joint distribution of the relative heights for the two boundary conditions. The form of these 

joint distributions will be used later in Section IV for the calculation of the MRH distribution in the stationary state. 



L 



A. Periodic Boundary Conditions 

We first consider the periodic boundary condition, H(Q,t) — H(L,t). The relative height h(x,t) — H(x,t) — 
J H(x,t)dx/L satisfies the same boundary condition, h(0,t) = h(L,t) and also satisfies the same evolution equation 
(28) as the actual height H(x,t). Since the relative height h(x,t) is a periodic function, one can decompose it into 
a Fourier series, h(x,t) = YZ=-ooK m ^) e2mmx,L - Substituting this in Eq. (28), one finds that different nonzero 
Fourier modes decouple from each other and one can easily calculate any correlation function. For example, the 
height-height correlation function is given by the exact expression, 

C(x = \ Xl -X 2 \,t,L) = (h( Xl ,t)h(X2,t)) = J-J2^- \l - e -^ m H/L^ e 2nim( Xl - X2 )/L^ (3Q) 

47r z m z L J 

where the sum ranges from m = — oo to m = oo excluding the m = term. 

In the growing regime when t « L 2 1 we can make a change of variable k — 2irm^/t/L. Since t << L 2 , dk — 2n>/t/L 
is small and hence k can be considered as a continuous variable. One can then replace the sum in eq. (30) by an 
integral over the k space. This integral can be easily performed and one gets, 
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C(x,t) = VtG(x/Vtj , 



(31) 



with the scaling function, G(y) = y §e y / 8 — | yerfc (j//-\/8), where erfc(a;) = e " is the complementary 

Error function. The scaling function G(y) starts at G(0) = \Z"2/tt and decays for large y as G(y) « ^^ 2 e~ y2 ^ 8 . The 

onsite variance grows as , (^ 2 (0)) = C(0,i) ~ i 1 / 2 indicating that the typical relative height in this regime grows as a 
power law, h <~ where the growth exponent /3 = 1/4 for the (1 + l)-dimcnsional EW interface. The height-height 
correlation between two sites, at a finite time t << L 2 and separated by a large distance 1 << x « L, decays as 
C(x,t) w 4 ^|^ 2 e~ x2 / u . Thus, in this growing regime, the correlations decay over a length scale £(i) <~ t 1 / 2 « L 
which is much smaller than the system size L. 

In the opposite stationary regime when t >> L 2 , one can drop the exponential factor in Eq. (30). The correlation 
function becomes time independent. Summing the resulting series in Eq. (30) one gets, 



C(x,L) = 



L 
12 



1 



6x 
T 



(32) 



Thus, the onsite variance (h 2 (0)) — C(0,L) = L/12. It is evident from Eq. (32) that in the stationary regime, the 
heights at different space points are strongly correlated. In Section V, we will calculate exactly the distribution of the 
maximum of these heights in the stationary regime. Thus, our result will provide an exact solution of the distribution 
of the extremum of a set of strongly correlated random variables. 

We are now ready to write down the joint distribution of heights P [{h}] in the stationary state. For simplicity, let 
us first consider the single site stationary height distribution P(h,t — > oo). Because of the linearity of Eq. (28) the 
stationary measure on h is Gaussian. Using (h 2 (0)) = C(0,L) = L/12, the stationary single site height distribution 
is Gaussian, 



P s t(h) = 




-6h 2 /L 



(33) 



Moreover, from Eq. (30), one can easily show that the slope-slope correlation function, (d x hd x >h) — > S(x-x') — 1/L'm 
the stationary state. The local slopes d x h are thus uncorrelated in the stationary state, except for the overall constraint 
due to the periodic boundary condition, L L dxd x h — 0, that gives rise to the residual l/L term. Collecting these 
facts together and remembering the constraint in Eq. (29), one readily writes down the stationary joint distribution 
of heights (a multivariate Gaussian distribution) , 



P [{h}} = A L e 



■i J L dr(d T hf 



5 [h(0) - h(L)} S 



f 

Jo 



h(r)di 



(34) 



where is a normalization constant and the two delta functions take care respectively of the periodic boundary 
condition h(0) = h(L) and the zero area constraint in Eq. (29). Note that in Eq. (34), r refers to the space (and not 
time) and varies over the interval < r < L. In fact, from now on we will use r to denote the space in order to make 
correspondence with the results derived in Sec. II. 

Before proceeding to calculate the normalization constant Al, we make one important remark here. The exponential 
factor inside the probability measure in Eq. (34) indicates that the stationary paths are locally Brownian, i.e., evolve 
in space as, dh(r)/dT — £(r), where £(r) is a Gaussian white noise with zero mean and a correlator, (£(t)£(t')) = 
S(t — t'). For the periodic boundary condition, the stationary path in space is, in fact, a Brownian bridge. This 
stationary Brownian measure also follows rather obviously even from the basic EW equation in Eq. (28) whose 
equilibrium Gibbs-Bolzmann state precisely corresponds to the exponential factor in Eq. (34). However, one has 



Jo K T ) dl 



in Eq. 



to be a bit more careful. The multiplicative factor representing the zero mode constraint, 5 
(34) indicates that only those Brownian bridges should be sampled which enclose a total zero area over the interval 
[0, L\. Normally, in writing down the stationary measure for the relative heights, one ignores this multiplicative delta 
function factor simply due to the fact that in most of the calculations on the interfaces, such as in the calculation of 
the distribution of the stationary width [30], this constraint does not play any important role. However, as we will 
see later in Section IV, this constraint does indeed play a major role in determining the maximal height distribution 
in the stationary state. Hence, it is worth being careful in writing the full exact stationary joint distribution of the 
relative heights keeping all the factors explicitly as done in Eq. (34). 
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The normalization constant Al : To determine the normalization constant, we note that if one integrates 
over the heights at all the intermediate points in Eq. (34) but keeping the heights at the two ends fixed at, say, 
h(0) = h(L) = u, one should recover the stationary single point height distribution P s t(^o = u ) m Eq. (33), .i.e., 



f hiL ^V h (r)e-iL L ^S [ 

Jh(0)=u JO 



h{r)di 




(35) 



To evaluate the path integral on the left hand side of Eq. (35), we use a simple and elegant method by identifying 
the path integral as the propagator of a random acceleration process that has been well studied [29]. To make the 
connection between the two problems, we first note that in the stationary state, the height h(r) locally evolves in 
space as a Brownian motion, dh(r)/dT — £(t) where £(r) is a Gaussian white noise with zero mean and a correlator, 
(£(t)£(t')) = 5(t — t'). But, as mentioned earlier, this by itself is not enough since it does not take into account the 



important delta function constraint 5 



Jo h (T)di 



To incorporate this special condition, let us define a new variable, 



X(t) = J r h{r)dr. The variable X(t) then evolves as, d 2 X/dr 2 = dh/dr = £(r) and hence can be identified as the 
position of a particle which is randomly accelerated. Thus, to take into account the full probability measure in Eq. 
(34), one has to consider the joint evolution, in space, of both the variables [X(r),/i(r)] where X(t) is the position 
and the h(r) = dX/dr is the velocity of the random accelerator. The joint propagator of this random acceleration 
process G[X, h, T\X , ho, 0], i.e., the probability that the process reaches its final value [X, h] in time T starting from 
its initial value [X , h ] at time 0, can be easily calculated and is well known [29] 



G[X,h,T\X ,ho,0] 



ttT 2 



exp 



~ (X - X - h T) (X-X - hT) - - h ) 2 



(36) 



The next step is to note that the path integral on the left hand side of Eq. (35) is just the propagator of the random 
acceleration problem in Eq. (36), going from its initial value [X n — 0, ho — u] to its final value [X = 0, h = u] in 
'time' L and hence is given by, 



r h ( L )= u i ri- , ,«. ,,2 T r L 

/ Vh{r)e—-h dT{a ^ 5 / h(r)dT 

Jh(0)=u JO 



G[0, u, L\0, u, 0] = -^-2 exp [-6u 2 /L] . 



(37) 



Substituting this result in Eq. (35) and cancelling out the exponential factors from both sides, one gets 

A L = V2^L 3 / 2 . (38) 



B. Free Boundary Conditions 

In this subsection, we consider the free boundary condition where the slopes of the interface vanish at the two ends, 
d x h = at x = and x — L. This boundary condition arises naturally if one considers a spatially discretized version 
of the continuum EW equation, 

dh ^^ = h(i + l,t) + h(i - 1, t) - 2h(i, t) + r)(i, t), (39) 
at 

where i = 1, 2 . . . , L and h(i,t) can be interpreted as the displacement of the i-th bead of a polymer chain where the 
beads are connected via harmonic springs, such as in the Rouse model [31]. If the chain forms a cycle, one has the 
periodic boundary condition. On the other hand, if the two ends of the chain are free, the beads at the two end points 
feel only one sided interactions, i.e., 

dk{L ' t] --h(L-l,t)-h(L,t) + r,(L,t), 



dt 
dh{l,t) 
Jt 



= h(2,t)-h(l,t) + V (0,t). (40) 



The end point equations (40) will be consistent with Eq. (39) provided one incorporates the boundary conditions, 
h(L,t) — h(L + l,t) and h(l,t) = h(0,t). In the continuum space, this boundary condition is equivalent to having 
d T h = at x — and x = L. 
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The zero slope conditions at the two end points indicate that one should decompose the height variable into a 
cosine series, h(x,t) — J2m=i h{m, t) cos(mirx/L). As in the periodic case, the different modes evolve independently 
and one can easily calculate the height-height correlation function. We find, 



2L 1 

C(x 1 ,x 2 ,t,L) = (h(x 1 ,t)h(x 2 ,t)) = — V* — 

7r * — ' m 

m—l 



^ _ e -2rn 2 Tr 2 t/L 2 



(■m-KX\\ (nvKX 2 \ 



(41) 



Note that, unlike the periodic case in Eq. (30), the correlation function C{x\, X2, t, L) is not translationally invariant 
and depends on both co-ordinates X\ and x 2 . 

In the growing regime when t << L 2 , one can again replace the sum in Eq. (41) by an integral over the k = irm/L 
space. The resulting integral can be easily performed and one finds, 



C(xi,X2, t) — V 



q [ \xi - x 2 \ \ +G f\xi +X2\ 



V 



V 



(42) 



with the same scaling function, G(y) = y §e y / 8 — \ ycrfc (y/V), as in the periodic case. Note that in the limit 
when both x\ and x 2 are large (but both small compared to L) with the distance between them x = \x\ — x 2 \ fixed, 
Eq. (42) reduces to the same correlation function as in the periodic case in Eq. (31). This is physically expected, 
since in the growing regime t << L 2 , the system is insensitive to the boundary conditions. 

In the stationary regime when t >> L 2 , the time independent correlation function can be obtained by dropping the 
exponential factor in Eq. (41) and summing the resulting series. We get, 



C(xi,x 2 ,L) = L 



\~^l^ Xl+ X2 + \ Xl ~ X2 ^ + 2I2 ( x i + X2 ) 



(43) 



In particular, the onsite variance, (h 2 (x)} — C(x,x,L) = L [| — (l — ^)] now depends on the site co-ordinate x 
since the translational invariance is lost with free boundary conditions. Note that, as in the periodic case, the heights 
are strongly correlated in the stationary state. 

The stationary single point height distribution is Gaussian, P st (h(x)) = e~ h ( X )/ 2C ( X , X > L ) j y/2nC(x, x, L). In par- 
ticular, the stationary height distribution at x = is given by 



(Mo)) 



2irL 



-3h 2 (0)/2L 



(44) 



One can also calculate the slope-slope correlation function from Eq. (41) and one finds that for t >> L 2 , (d x hd x >h) — > 
5(x — x'), implying that the slopes get completely uncorrelated in the stationary state. Once again, collecting these 
facts together and taking into account the constraint in Eq. (29), one finds that the stationary joint distribution of 
heights for the free boundary case, for fixed boundary heights h(0) — u and h(L) — v, is given by 



P[{h} lU ,v] = B L e- l2 Io dT ^ h)2 



h(r)dT 



S[h(0)-u] 6[h(L)-v] 



(45) 



where the normalization constant B L is yet to be determined. 

The normalization constant B L : The constant B L is calculated using the similar method as in the periodic 
case. First, we integrate over the heights at all the intermediate points in Eq. (45) but keeping the heights at the 
two ends fixed, say h(0) — u and h(L) — v. This gives us the joint stationary distribution of heights at the two end 
points, P st [u,v]. If we now further integrate over v, we should get the single point height distribution P st [h(0) = u] 
given in Eq. (44). Thus we have, 



B r 



r h(L)=v 

dv / Vh(r)e 

Jh(0)=u 



i f o L dr{d T h) 2 



I' 

JO 



h{r)dT 



2irL 



-3u 2 /2L 



(46) 



As in the periodic case, we now identify the path integral on the left hand side of Eq. (46) as the propagator 
G[0,v,L\0, u, 0] of the random acceleration process in Eq. (36) going from the initial point [0, u] to the final point 
[0,v] in 'time' L. This gives, 
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ML)=V 

/ T>h{ 

Jh(0)=u 



r)e -i./>^) 2 



h(r)dT 



G[0,v,L|0,u,0] = Xr* exp 

71"L Z 



6 2 2 
-— uv — -(u — v) 
Li hi 



(47) 



We substitute this expression of the path integral on the left hand side of Eq. (46) and then integrate over v. 
Cancelling the exponential factors from both sides of Eq. (46) gives 



B L = L. 



(48) 



IV. MAXIMAL RELATIVE HEIGHT DISTRIBUTION IN THE STATIONARY REGIME: 

EDWARDS- WILKINSON INTERFACE 

In this section, we derive the principal new results of this paper, namely the exact probability distribution of the 
maximal relative height (MRU) in the stationary state of the Edwards- Wilkinson equation in a finite one dimensional 
system of size L, both for periodic and free boundary conditions. 



A. Periodic Boundary Conditions 

Consider the Edwards- Wilkinson model in Eq. (28) in the stationary regime when t >> L 2 . The exact form of 
the stationary joint distribution of the relative heights P [{h}] is given in Eq. (34) with the normalization constant, 
Al = \/2nL 3 / 2 . We wish to compute the probability distribution of the MRH, max [{ft,}]. Let us first define the 
cumulative distribution of the MRH, F(h m , L) = Prob [max{/i} < h m , L}. The distribution of the MRH is simply the 
derivative, P(h m ,L) = dF(h m , L)/dh m . Clearly F(h m ,L) is also the probability that the heights at all points in 
[0, L] are less than h m and hence can be written using the measure in Eq. (34) as a path integral, 



/fT-m /•h(L)=U 1 C L 2 f U 

du / m(r)e~ 5 Jo dT ^ h) S / h(r)dT 

-oo Jh(0)=u JO 



I(h m ,L), 



(49) 



where I(h m , L) = Y[ T =o @{h"m — M r )) i s an indicator function which is 1 if all the heights are less than h m and zero 
otherwise. Due to the periodic boundary condition, h(0) = h{L), all the paths inside the path integral propagate from 
their initial value h(0) = u to their final value h(L) = u. Note also that one needs to integrate over the boundary 
value h(0) — h(L) — u over the allowed range — oo < u < h m . The upper bound u < h m follows from the fact that, 
by definition, h m is the maximum over all heights in [0, L] including the two end points. 

The next crucial step is to make a change of variables, x(t) = h m — h(r) and u' = h m — u in the path integral in 
Eq. (49) which gives, 



POO px(L) — U 1 P L 2 r 1 ^ 

F(h m ,L) = A L du' Vx(T)e~*J« dT(dTX) S / x(r)dT - A 

Jo Jx(0)=u' Jo 



I(L), 



(50) 



where I(L) = J^ T=0 9(x(t)) and A — h m L. Physically, this change of variable just means that we are now measuring 
the interface height relative to the maximum and denoting it by —x(t). The maximal height now corresponds to the 
level x = 0. Note that h m appears only through the quantity A in the delta function, and hence F(h m , L) = J 7 (A, L). 
In the subsequent analysis, we will keep a general A in Eq. (50) and will substitute A = h m L only in the final formula. 

We already start seeing the formal similarity between Eq. (50) and Eq. (8) in Section II for the distribution of area 
under a Brownian excursion. To be more precise, we take the Laplace transform with respect to A in Eq. (50) which 
gives, 



f 

Jo 



T{A,L)e- XA dA = A L 



Jo 



du' < u'\e- HlL \u' 



>= A L Tr 



-H X L 



(51) 



where Tr is the trace and H\ is precisely the same Hamiltonian that appeared in Section II namely, H\ = — ^^jt+^i^) 
where V\(x) is a triangular potential: V\(x) = Xx for x > and V\{x) — oo for x < 0. The Hamiltonian Hi has 
only discrete eigenvalues that were obtained in Section II, E k = a fe A 2 / 3 2~ 1//3 for k — 1,2, .. ., where ct^'s are the 
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magnitude of the zeros of Ai(z) on the negative real axis. The trace in Eq. (51) can be easily evaluated knowing 
these eigenvalues and one gets using Al — V2n L 3 ' 2 , 



poo ^ 

/ T(A,L)e- XA dA = V2^L 3 / 2 ]Te 
J ° k=i 



OO 

akX 2/3 2 -l/3 L 



(52) 



Next we invert the Laplace transform in Eq. (52) using Bromwitch formula and substituting A = h m L we get the 
cumulative MRH distribution, 

F(h m ,L) = V2^L 3 / 2 ^ e Xh mL y e - ak x^2-^L^ (53) 

where the integration is along any imaginary axis whose real part Ao must be to the right of all singularities of the 
integrand. Taking derivative with respect to h m in Eq. (53) and making a change of variable, A = sL~ 3 / 2 , we arrive 
at our main result, P(h m , L) = L~ x / 2 f {h m L~ x l 2 ) for all L, where the Laplace transform of f(x) is given by 



/• oo 

Jo k=i 



-^ 2/32 - 1/3 . (54) 



Comparing with Eq. (13) in Section II, we identify the scaling function f(x) as the Airy distribution function that 
characterizes the probability distribution of the area a under a Brownian excursion over a unit interval. 

Thus, the main result of this subsection has been to establish an equivalence in law between two rather different 
objects: (i) the scaled maximal relative height h m j\TL in the stationary state of the EW equation (28) over a linear 
substrate of size L and (ii) the area a under a Brownian excursion over the unit time interval [0, 1], 

tI 55 fl ' (55) 

where = means that the object on the left hand side has the same probability distribution as the object on its right 
hand side. 

Moments of the MRH: Using the scaling form, P(h m ,L) = L' 1 / 2 f (hmL- 1 / 2 ) for all L, we get (h m ) = M n L n ' 2 
where M n = f(x)x n dx are the moments of the Airy distribution function that were computed in Eq. (16) of 
Section II. In the context of the MRH, only the second moment {h m ) = 5L/12 was computed previously in [21] by 
using a different method. Here, we can compute all the moments of MRH recursively. 

Asymptotic tails of the MRH distribution: In order to determine the asymptotic behaviors of the MRH 
distribution, P(h m ,L) = L -1 / 2 / (hmL^ 1 / 2 ) , we need to know the tails of the function f(x). As mentioned in the 
introduction, Takacs was able to formally invert the Laplace transform in Eq. (1) and obtained an expression of f(x) 
in terms of confluent hypergeometric functions as given in Eq. (2). It is easy to obtain the small x behavior of f(x) 
from Eq. (2), since only the k = 1 term dominates the sum near x — > 0. Using U(a, b, z) ~ z~ a for large z [4], we get 
as x — > 0, 



tt \ 5 9 / 2 -5 

f(x) -» ^j-<V x exp 



2a 



27x 2 



(56) 



where ol\ = 2.3381074 .... In the context of the MRH distribution, this essential singular tail of the scaling function 
near x — > was conjectured in [21] based on numerics, though the exact form was missing. The asymptotic behavior 
at large x is more tricky to derive [5] from Eq. (2). However, it is possible to guess the leading behavior of f(x) for 
large x by examining the behavior of moments M n . Using the recursion relations for the moments M n in Section 
II (Eq. (15)), one finds that for large n, M n <~ (n/12e)"/ 2 . Substituting an anticipated large x decay of the form, 
f(x) ~ cxp[— ax b ] in M n — J °° f(x)x n dx, we get M n <~ (n/abe) n / h for large n. Comparing this with the exact large 
n behavior of M n we get a = 6 and 6 = 2, indicating a Gaussian tail as x — > oo, 

f(x)~e- 6 * 2 . (57) 
However, we are not able to determine the amplitude (or the next subleading correction) of this Gaussian tail. 
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FIG. 4. The scaling function f(x) associated with the MRH distribution , P(h m , L) = Z,~ 1//2 f(h m /y/L) for the EW equation 
with periodic boundary condition. The numerical curves (shown by symbols) are obtained by collapsing the data obtained from 
the numerical integration of Eq. (58) for three system sizes L — 256 (circles), L — 384 (squares), and L = 512 (diamonds). 
Also plotted is the Mathematica generated analytical scaling function in Eq. (2) as shown by the solid line. 

Comparison with numerical simulations: Finally, we would like to compare our exact result for the MRH 
distribution, P(h m ,L) = L" 1 / 2 / (h m L- 1/2 ) with f(x) given in Eq. (2), to the results of numerical simulations. We 
numerically integrated the space-time discretized form of Eq. (28) , 

H(i, t + At)- H(i, t) = At [H(i + 1, t) + H(i - 1, t) - 2H(i, t)\ + 7j t (t)V2At, (58) 

where r)i(t)'s are independent and identically distributed random variables for each i and t and each drawn from a 
Gaussian distribution with zero mean and unit variance. We used periodic boundary condition, i.e., -ff (0, t) = H(L, t) 
and H{L + l,t) = H(l,t). We chose the time step At = 0.01, and we checked that the results were stable. We used 
different system sizes L = 256, L — 384, and L — 512. We first evolved Eq. (58) for a long time (typically 2 x 10 6 
Monte Carlo time steps) and ensured that the system has reached the stationary state. After that, the system was 
further evolved and the data for the histogram of the MRH was sampled at times which are far apart from each other 
to avoid correlations between them. The number of samples used were typically 2 x 10 6 . In Fig. 4, we have plotted 
the numerical scaling function generated by collapsing the data for the three system sizes and also plotted the Airy 
distribution function f(x) by evaluating the sum in Eq. (2) using the Mathematica. The agreement between the 
analytical and the numerical scaling function is very good. 



B. Free Boundary Conditions 



In this subsection we compute the stationary MRH distribution of the EW equation (28) with free boundary 
conditions, d x h — at x — and x = L. The stationary joint distribution of relative heights P [{h},u,v] (with 
boundary heights fixed at h(0) — u and h(L) — v) is given in Eq. (45) with Bl = L. As in the periodic case, we 
define the cumulative MRH distribution, F(h m ,L) = Prob [max{/i} < h m ,L}. Since F(h m ,L) is also the probability 
that the heights at all points in [0, L] are less than h m , one can use Eq. (45) to write 



F(h m ,L) = L 



du 



dv 



h(L)=v 



h(0)=u 



h(r)dT 



I(h m ,L) 



(59) 



where I(h m ,L) = Y[t=o^(^"i — M t ))j as before, is an indicator function which is 1 if all the heights are less than 
h m and zero otherwise. Note that the boundary heights u and v are integrated over their common allowed range 
[-oo, h m }. 
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The next step is to make the change of variables, x(r) = h m — h(r), u' = h rn — u and v' — h rn — v where —x(t) 
now denotes the interface height measured relative to the maximal height. The maximal height now corresponds to 
the level x = 0. With these change of variables, Eq. (59) now reads 



F(h m , L) = L du' dv' Vx(T)e~iJ° dT{dTX > 5 / x( T )d T 

JO JO Jx(0)=u' Jo 



I(L) 



(60) 



where I(L) = Y\ T=0 9(x(t)) and A = h m L. The only dependence on h 7n in Eq. (60) appears through the quantity 



A = h m L. Hence, F(h m , L) = J 7 (A, L). Taking Laplace transform with respect to A in Eq. (60) gives, 

f-OO f-OO /'OO 

/ F(A,L)e- XA dA = L du' dv' < u'\e~^ L \v' >, 
Jo Jo Jo 



(61) 



where the Hamiltonian Hi is the same as in the periodic case. Eq. (61) for the free boundary case should be 
compared to the corresponding result for the periodic case in Eq. (51) where the integral over the common end 
points became a trace due to the periodic boundary condition. In the present case, we do not have a trace and 
evaluating the integral in Eq. (61) would require the knowledge of both the eigenfunctions and the eigenvalues of H\. 
Fortunately, the eigenfunctions can be determined exactly in terms of the Airy function as explained in Section II. 
Expanding the propagator < u'\e~ HlL \v' > in Eq. (61) in the energy eigenbasis of Hi and using the exact eigenvalues 
Ek = afeA 2 / 3 2 -1 / 3 for k = 1, 2, . . . (afe's being the magnitude of the zeros of Ai(z) on the negative real axis) we get, 



I"* F{A,L)t 
Jo 



-XA 



dA 



k=i 



p OO pOO 

Jo Jo 



Substituting the exact form of the eigenfunctions Vfe (x) from Eq. (9) and simplifying, we get 

j- OO 



I 



HA, L)e~ XA dA = — !— CK) * 
[ZA > k=i 



(62) 



(63) 



The constant C{ak) is given by the expression, 



C{a k ) 



J x a Ai(z)dz J™ a Ai{z)dz 



r^(z)dz 



-oc k 

pOC 



Ai>i(-a k ) 



(64) 



where we have again used the identity [28], J^° a Ai 2 (z)dz = [Ai'(— afc)] 2 . Note that C(ak) = B 2 (otk), where B(ctk) 
is given in Eq. (21). The constants C(ak)'s can be determined numerically very precisely from Eq. (64) using the 
Mathcmatica. As examples, we have determined the first few values in Table I. Note that C(afc)'s oscillate as k 
increases. 



k 




C{a k ) 


1 


2.33810... 


3.30279... 


2 


4.08794 . . . 


1.01282... 


3 


5.52055 . . . 


1.78218... 


4 


6.78670 . . . 


0.90519... 


5 


7.94413 . . . 


1.39473... 


6 


9.02265 . . . 


0.83181... 


7 


10.04017... 


1.19915... 


8 


11.00852... 


0.77846 . . . 


9 


11.93601... 


1.07584... 


10 


12.82877... 


0.73725 . . . 



Table 1. The second column shows the magnitude at of the fc-th zero of the Airy function Ai(z) on the negative real 
axis upto k — 10, obtained using the Mathematica. The third column shows the corresponding values of the constant 
C(afe) in Eq. (64), again generated via the Mathematica. 

We formally invert the Laplace transform in Eq. (63), take the derivative with respect to h m and then make a 
change of variable, A = sL -3 / 2 . This leads us to the exact MRH distribution for the free boundary case, P(h m , L) = 
£-1/2 p ^ m £-1/2^ f or a n ^ where the Laplace transform of the scaling function F(x) is given by, 
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P oo 00 

/ F{x)e- sx dx^2- 1 ^s 2 / 3 J2 C ( a k)' 
Ja k=i 



(65) 



with C(afe)'s determined from Eq. (64). One can check from Eq. (65) that the function F(x) is normalized, 
f °°F(x)dx = l. 

In order to make contact with the area under a Brownian meander discussed in Section II-B, we define a new 
variable t = 2 _1 / 3 s 2/3 and rewrite Eq. (65) as, 



/>oo 

/ F(x 
Jo 



~dx = E 



OO 



(66) 



where E denotes the expectation value with respect to the distribution F(x). Taking a further Laplace transform of 
Eq. (66) with respect to t we get 



P oo 

/ e~ ut E 




dt 


'0 







C{a k ) _[CMz)dz] 



—J (u + a k f 



Ai(u) 



(67) 



The second equality in Eq. (67) is proved in Appendix-A. Comparing Eq. (67) with Eq. (24) of Section II-B, one 
finds that the former is just the square of the latter. Using the result in Eq. (24) and the convolution theorem of 
Laplace transform, we can then invert the Laplace transform in Eq. (67) to obtain the following result, 



E 




-L 









dt' 



/O TTy/t'(t - f) 



E 



-V2*' 3/2 ai 



E 



-V2(t-t') 3/2 a 2 



(68) 



where a\ and 02 are the areas under two independent Brownian meanders each over unit intervals. The result in Eq. 
(68) is valid for all t. In particular, putting t = L and x = h m /\/L, we get 



E 



-V2Lh„ 



dt' 



ny/t'(L - t') 



:E 



E 



-V2(L-t') 3/2 Q2 



(69) 



which has a very nice and simple physical interpretation. Consider the stationary interface profile over the interval 
[0,L] after we have made the shift with respect to the maximum, x(t) = h m — h(r). A typical picture is shown in 
Fig. 5. 




Of L 
FIG. 5. A stationary profile of the interface x(t) — h m — h(r) with the free boundary condition as seen from the maximum 
h m located at t' . The point t' divides the path into two statistically independent Brownian meanders with areas Ai and A 2 
respectively. 

Let the location of the maximum be denoted by t'. Note that the location < t' < L of the maximum of a free 
Brownian motion over the interval [0, L] is a random variable whose normalized probability distribution P(i') is well 
known [32], 



P(f) 



ny/t'(L - t') 



(70) 



which also describes, in the context of the interface with free boundary conditions, the distribution of the position of 
the maximum of the interface. It is clear from Fig. 5 that the maximum h m at t' (or rather the minimum x = in Fig. 
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5) divides the path into two Brownian meanders, one over the left interval [0, t'] and the other over the right interval 
[t', L] whose interval length is (L — t'). Since a Brownian motion is Markovian, these two meanders are statistically 
independent of each other. Thus, the areas under the two meanders, A\ — ait' 3 / 2 and A 2 = a 2 (L — t') z / 2 are also 
statistically independent, except that their joining point t' is distributed according to Eq. (70). Thus, Eq. (69) has 
now a probabilistic interpretation, 



hm.L = 



[A 1 (t') + A 2 {L-t')}P(t')dt', 



(71) 



where P(t') is given in Eq. (70). The random variable h m L has therefore the same law as the sum of the two meander 
areas ( with the location of their joining point t' integrated over with the measure P{t') in Eq. (70)). This result in 
Eq. (71) for the free boundary condition is thus more involved than its periodic counterpart in Eq. (55). 

Moments of the MRH: Using the scaling form, P(h m ,L) = L" 1 / 2 F(h m Z," 1 / 2 ) for all L, we get (/i™ ) = ^L"/ 2 , 
where fi n — J °° F{x)x n dx are the moments of the function F(x) defined in Eq. (65). The computation of yu„'s from 
Eq. (65) is nontrivial. To make progress, we make use of Eq. (68). We note that the right hand side of Eq. (68) 
involve the areas cti and a 2 of two independent Brownian meanders over unit intervals. Now, the moments associated 
with the area of a Brownian meander over a unit interval a n = E[a n ] were computed by a recursive procedure in 
Section II-B and the results are given in Eq. (27). By expanding the exponentials on the right hand side of Eq. (68) 
in power series and then performing the integral over t' we get 



E 



(-V2)" 



^ oc oc 

* r ( n l + : ) r ( n2 + 

ni=0ri2=0 



^,321 + 1^3^+1^ anian2t s {ni+n2)/ ^ 



(72) 



where B(x,y) is the standard Beta function and the moments a n 's of the area under a single meander over unit 
interval were computed in Eq. (27) in Section II-B. Next, expanding the left hand side of Eq. (72) in a power series 
and matching identical powers of t, we can express the moment /i„ = J °° F(x)x n dx in terms of the a„'s, 



1 " 
= - V 



7%\ ( 3m + 1 3(n — m) + 1 
TT ^ n \ m j \ 2 2 

m— 



(73) 



Using the known values of a„'s from Eq. (27), one can then recursively determine the moments fi n from Eq. (73). 
For example, the first few values are 



Mo = 1, 



Mi = 



M2 = 



17 

24' 



M3 = 



123 

140 V TT 



(74) 



Asymptotic tails of the MRH distribution: To determine the asymptotic tails of the MRH distribution, 
P(h m , L) = L~ X / 2 F (h rn L~ x / 2 }, we need to know the asymptotic behavior of the scaling function F(x) defined in Eq. 
(65) in the limits x — > and x — > oo. Following similar method as Takacs in the periodic case in Eq. (2), we were able 
to invert the Laplace transform in Eq. (65) and express it in terms of the hypergeometric function. The derivation is 
presented in Appendix-B. We get, 



F(x) 



2-1/3 

3^x 7 / 3 



J2 a k C(a k ) [U (1/6, 4/3, b k /x 2 ) + 2U (-5/6, 4/3, b k /x 2 )] e 



(75) 



fe=i 



where b k = 2a|/27, C(a k ) is defined in Eq. (64) and U(a,b,z) is the confluent hypergeometric function. Thus, for 
the free boundary condition, the scaling function F(x) is different from that of its periodic counterpart namely the 
Airy distribution function in Eq. (2). The function F(x) appears to be new and hence it lacks a name. We call this 
function F(x) as the i^-Airy distribution function, where F refers to the free boundary condition. 
In the limit x — > 0, only the k = 1 term dominates the sum in Eq. (75) and we get 



F(x) 



2V2 
27V^ 



\ 7/2 -4 
C(ai)a 1 ' x exp 



2a? 
~27x 2 



(76) 



where ct\ — 2.33810 . . . and C{a\) — 3.30278 . . . from Table I, evaluated using the Mathematica. Thus the function 
F(x) decays slightly faster as x — > than its periodic counterpart in Eq. (56). Finding the asymptotic behavior 
of F(x) as x — > 00 turns out to be more tricky. However, as in the periodic case, the leading large x behavior can 
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be guessed by analyzing the moment \i n for large n. It is not difficult to analyze Eq. (73) for large n, knowing the 
properties of a n 's. We find that for large n, fi n ~ [n/2e] n ^ 2 . Following exactly the same procedure as in the periodic 
case, we find that F(x) has a leading Gaussian tail for large x, 



F(x) 



-3x 2 /2 



(77) 



that falls off less rapidly than the periodic case f(x) ~ e~ 6x in Eq. (57). 

Comparison with numerical simulations: Our exact analytical result for the MRH distribution for the free 
boundary condition is: P(h m ,L) = L~ X I 2 F (/i m L -1 / 2 ) where the scaling function F(x) is given in Eq. (75). We 
have evaluated the sum in Eq. (75) using the Mathematica. It turns out that the sum is rapidly convergent and it is 
sufficient to keep terms upto k = 10 in the sum. We also numerically integrated the discretized Edwards- Wilkinson 
Eq. (58) with free boundary condition, i.e., H(0,t) — H(l,t) and H(L + l,t) = H(L,t) and calculated the MRH 
distribution. In Fig. 6, we compare the analytical scaling function to the one obtained by numerical integration. The 
number of steps used to reach the stationary state as well as the number of sample averages were the same as in 
the case of the periodic boundary condition. The numerical data shown in Fig. (6) were obtained by collapsing the 
histograms for three system sizes L = 256, L — 384, and L — 512. The agreement is again very good. 




FIG. 6. The scaling function F(x) associated with the MRH distribution , P(h m , L) = L" 1 ^ 2 F(h m /VL) for the EW equation 
with free boundary condition. The numerical curves (shown by symbols) are obtained by collapsing the data obtained from the 
numerical integration of Eq. (58) for three system sizes L = 256 (circles), L — 384 (squares), and L — 512 (diamonds). Also 
plotted is the Mathematica generated analytical scaling function in Eq. (75) as shown by the solid line. 

V. MAXIMAL RELATIVE HEIGHT DISTRIBUTION IN THE STATIONARY REGIME: 

KARDAR-PARISI-ZHANG INTERFACE 



In the previous section, we derived the exact distribution of the MRH in the stationary regime of the (1 + 1)- 
dimensional Edwards- Wilkinson interface evolving via Eq. (28). A natural next step would be to extend our method 
to other (1 + l)-dimensional interfaces. In this section, we study the MRH distribution in the stationary regime of 
another very well studied interface that evolves via the nonlinear KPZ equation [23], 



dH(x, t) _ d 2 H(x, t) ( dH(x, t) 
+ AI 



dt 



dx 2 



dx 



ri(x,t), 



(78) 



where rj(x,t) is a Gaussian white noise with zero mean and a correlator, (r](x,t)rj(x't')) — 26 (x — x')5(t — t'). As 
before, we subtract the zero mode and focus on the relative height, h{x,t) — H(x,t) — H(x,t) whose distribution 
reaches a stationary state in the long time limit in a finite system. 

It is well known that the dynamical exponent is z = 3/2 for the (1 + l)-dimensional KPZ equation [16]. Thus, 
we expect that for times t » L 3 / 2 , the system will reach a stationary state. We wish to compute the distribution 
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P(h m ,L) of the MRH, h m = max [{ft,}], in this stationary state. In order to compute this distribution, we need to 
know the stationary joint distribution of the relative heights P [{h}] for the KPZ equation. In an infinite system, it is 
well known [16] that in (1 + l)-dimcnsions, the joint distribution of heights approaches a stationary state as t — > oo, 
where P [{H}\ oc exp [— 5 / {d x H) 2 dx\ . This is a somewhat surprising result since the dependence on the nonlinear 
term completely drops out of the stationary measure. Indeed, this turns out to be a very special property only valid 
in (1 + l)-dimcnsions. This result can be proved by writing down the full Fokker-Planck equation for the joint height 
distribution P [{H}, t] and then showing explicitly that indeed P [{H}] oc exp \—\ J (d x H) 2 dx\ is a stationary solution 
of the Fokker-Planck equation [16]. In proving the last step, one neglects various boundary terms that arise out of 
integration by parts [16] which is justified only in the limit of an infinite system. However, in a finite system, the 
boundary terms are usually nonzero and it is no longer easy to find the stationary measure. 

The only exception seems to be the periodic boundary case where, even in a finite system, the boundary terms can 
be shown to vanish. Hence in this case the joint distribution of the relative heights P [{h}] for the KPZ equation has 
the same expression as in the EW case, and is given by Eq. (37) with the normalization constant Al = \/~2tt L 3 / 2 . 
Therefore, for the periodic boundary condition, we expect that the stationary MRH distribution for the KPZ equation 
will also be identical to that of the EW case, i.e., P(h m ,L) = L~ x / 2 f (/i m L~ 1,/2 ) where the scaling function f(x) is 
the Airy distribution function in Eq. (2). 

An interesting challenge is to verify this analytical prediction numerically by directly integrating the KPZ equation 
(78). The problem is to find an appropriate spatial discretization scheme that will correctly represent the nonlinear 
term in the continuum equation (78). Most of the past studies have used the following natural choice or its simple 
variants [33-35], 



H(i,t + At) - H(i,t) = At 



(H (i + 1, t) + H (i - 1, t) - 2H (i, t)) + ^(H(i + l,t)-H(i- 1, t)f 




where 77^ (£) 's are independent and identically distributed random variables for each i and t and each drawn from a 
Gaussian distribution with zero mean and unit variance. However, it is well known [35,36] that certain properties 
of Eq. (79) are rather unphysical and are fundamentally different from those of the continum equation (78). In 
other words, the discrete equation (79) fails to capture the correct continuum KPZ evolution [37]. In this paper, we 
have instead used a discretization scheme proposed by Lam and Shin [38] which circumvents the problems mentioned 
above. According to this scheme, one uses the equation (79) except that the nonlinear term in Eq. (79) is replaced 
by the following expression [38], 



(H {i + l,t)-H (i, t)f + {H{i + l,t)-H (i, t)) (H (*, t)-H(i-l, *)) + (H (*, t) - H (i - 1, t)f . (80) 



The advantage with this scheme is that one can prove analytically that the Fokker-Planck equation associ- 
ated with this discrete model admits a stationary solution for the periodic boundary case [38], P [{H}] oc 

exp — i J2i=i (H (* + 1) — H W) 2 which, in the contiuum limit, correctly reproduces the stationary measure of 

the continuum KPZ equation, i.e., P [{H}] oc exp — | (d T H) 2 dT . We note that a similar but slightly different 

scheme that also leads to the Gaussian stationary state has been proposed earlier [39] . 

We evolved the discretized equation (79) (but with the nonlinear term in Eq. (79) replaced by the Lam-Shin 
nonlinear term in Eq. (80)) with periodic boundary condition and chose A = 1 and At = 0.01. The steady state 
MRH distribution P(h m , L) was obtained for three different system sizes L = 256, L = 384 and L = 512. The three 
data sets were collapsed onto a single scaling plot, P(h m , L) = L~ x / 2 f (/i m L~ 1//2 ) . The numerically obtained scaling 
function is then compared in Fig. 7 with the analytical scaling function f(x) given in Eq. (2). The agreement is 
excellent. Thus, for the periodic boundary condition, the stationary MRH distribution of the (1 + l)-dimensional 
KPZ equation is also described by the Airy distribution function, as in the case of the EW equation. However, we 
expect that this superuniversality of the Airy distribution function holds only in (1 + l)-dimensions, and not in higher 
dimensions. 
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FIG. 7. The scaling function f(x) associated with the MRH distribution , P(h m ,L) — L -1 ' 2 f{h m /yL) for the KPZ equation 
with periodic boundary condition. The numerical curves (shown by symbols) are obtained by collapsing the data obtained 
from the numerical integration of Eq. (79) (but with the nonlinear term replaced by the Lam-Shin term from Eq. (80)) for 
three system sizes L — 256 (circles), L = 384 (squares), and L = 512 (diamonds). Also plotted is the Mathematica generated 
analytical scaling function in Eq. (2) as shown by the solid line. 

We now turn to the case of the free boundary condition, i.e., with Neumann boundary conditions d x h = at the 
two ends x = and x = L. In the periodic case P[{H}} = cxp [—4 J (d x H) 2 dx] is explicitly the stationary solution 
of the Fokker- Planck equation. However, for the free boundary case, the same measure is no longer a stationary 
solution of the Fokker-Planck equation. Thus, unlike in the periodic case where the Airy distribution describing the 
MRH distribution for the Edwards- Wilkinson equation also describes the MRH distribution of the KPZ equation, for 
the free boundary case there is no apriori reason to believe that the F-Airy distribution for the Edwards- Wilkinson 
case will carry over to the KPZ case as well. In fact, as far as we know, the stationary measure for the continuous 
KPZ equation with the Neumann boundary condition is still not known explicitly. Thus, calculating analytically 
the corresponding MRH distribution remains an interesting open problem. Let us also remark about the numerical 
simulation with free boundary condition. As discussed before in this section, to simulate the continuous KPZ equation 
one needs to use an appropriate discretization scheme that will reproduce the correct continuum stationary measure. 
Arbitrary discretization scheme can lead to spurious results. For the periodic boundary condition, such discretization 
schemes are available thanks to the fact that one knows the stationary state explicitly. For the free boundary case 
where one does not know the stationary state explicitly, it is apriori not clear what type of discretization schemes one 
should use. 

In fact, to study the effects of different boundary conditions on the MRH distribution in nonlinear KPZ type 
interfaces, it may be more appropriate and easier to study the discrete growth models that belong to the KPZ 
universality class, rather than the continuous KPZ equation itself. Several such models are known [40], an example 
being the single-step model [41]. This single-step model can further be mapped to the asymmetric exclusion process 
(ASEP). A particle at a site i in ASEP corresponds to a downward slope of the interface height on one unit and a hole 
in ASEP corresponds to an upward slope of the interface step of one unit. The height hi(t) at site i in the interface 
model is thus related to the occupancy Ti in ASEP via the simple relation, hi + i(t) — hi(t) = 1 — 2rj(t) where Tj = 1 
(or 0) if the site i is occupied (or empty) in the ASEP. Over the past few years there have been extensive studies 
of ASEP with open boundary conditions where a particle enters the lattice through its left end at rate a and leaves 
the lattice throught its right boundary at rate j3 [42] . Several steady state properties of this system are known [43] . 
The open boundary conditions in ASEP means special growth rules at the boundaries of the interface model. Thus, 
many results from the ASEP can then be used to predict the properties of the interface model with special boundary 
conditions [43]. In particular, for a < 1, (3 < 1 and on the line a + (3 = 1, ASEP is known to have a factorized steady 
state, i.e., the Tj's become independent from site to site, each with a bivariate distribution p(r) = aS Tt i + (36 T fi. This 
means that the slope variables in the corresponding interface model also get uncorrelated in the stationary state. The 
stationary profile in space is then described by a one dimension random walk with drift, hi+i — hi = & where £j's are 
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i.i.d variables which take the value 1 (with probability [3 — 1 — a) and (with probability a). Thus, we have a random 
walk with a drift p = 1 — 2a. For the zero drift case, a = (3 = 1/2, we would then expect that the F-Airy distribution 
function, upto a scale factor, would describe the MRH distribution of the interface in the continuum limit. It would 
be interesting to extend the method presented here to calculate the MRH distribution for p, ^ 0. 

VI. MAXIMAL RELATIVE HEIGHT DISTRIBUTION OF INTERFACES IN THE GROWING REGIME 

In this section, we discuss the MRH distribution of a (1 + l)-dimensional interface in its growing regime, i.e., when 
the time 1 << t << L z where z is the dynamical exponent associated with the growth. The MRH distribution in 
the growing regime turns out to be more universal than its stationary counterpart and holds for generic interfaces. 
Logically, one would have preferred to discuss the growing regime (early time) before the stationary regime (late time). 
However, we will need to use some of the results for the stationary regime derived in Section IV in our discussion of 
the MRH distribution in the growing regime. Hence is the reversal of the order. The MRH distribution in the growing 
regime has been discussed previously in Ref. [21] and we will be using a similar line of arguments. Nevertheless, 
we decided to include this section partly because our results would be more precise and partly to make this paper 
self-contained and complete. 

Let us consider a growing interface starting from a flat initial condition at t = 0. Initially the heights are completely 
uncorrelated. As time grows, the correlation between heights also grow. As discussed in Section III, at any finite time 
t, the height-height correlation function decays with spatial distance over a characteristic correlation length £ (t) which 
grows algebraically with time, £(t) ~ t x l z . In the growing regime when the time 1 << t « L z , this correlation length 
is much smaller than the system size, £(t) << L. Thus, one can break up the whole system of size L into N = L/£ 
blocks each of size as shown in Fig. 8. The heights (or rather the relative heights) at two points belonging to 
two different blocks are basically uncorrelated since the distance between them is bigger than the correlation length 
£(t). However, inside a given block, the heights are strongly correlated. Let us denote the maximal relative height in 
block i by h m (i). These block MRH's are shown by the black dots in Fig. 8. 




Ox L 

FIG. 8. An interface in the growing regime with correlation length The system is divided into N — L/£(t) blocks each 

of length £(t). The heights in different blocks are essentially uncorrelated. The maximal relative height in each block is denoted 
by a black dot. Also shown is the global maximal relative height h m . 

We are interested in the distribution of the global MRH, h m = max [h m (l), h m (2), . . . , h m (N)\ where N is the number 
of blocks. Since the block variables /i m (z)'s are uncorrelated, we are thus interested in the extreme value statistics 
of a set of N uncorrelated random variables, a subject that has been rather well studied [17]. Furthermore, since 
£(t) « L, the system does not feel the presence of the boundaries and there is a translational invariance over the block 
index i. This indicates that the random variables /i TO (i)'s associated with different blocks are not just independent, 
but share the same distribution as well. Let us denote this common distribution by p(x,£) — Prob[/i TO (i) = 
where £ is the size of the block. Due to the independence of different blocks, the joint distribution P [{h m (i)}] of the 
block variables /i m (z)'s factorize, P [{h m (i)}] = Y\i P {h m {i), £). 

To guess the form of the distribution p(x, £) = Prob [h m (i) = x, £], one makes a reasonable assumption that within 
a given block of size £, the heights have already reached the stationary state even though the full system is still in the 
growing regime. Then p(x) is just given by the stationary MRH distribution with a system size £ and one can then 
directly use the results obtained in the previous section, namely the scaling form, p(x,£) = ^~ X I 2 W (x£~ 1//2 ). More 
generally, one should use the scaling form, p(x,£) — £~ a W {x£~ a ), where a is the roughness exponent of the surface. 
For the (1 + l)-dimensional EW and the KPZ equation, one has a = 1/2. The roughness exponent a is related to 
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the dynamical exponent z and the growth exponent (3 via the simple scaling law, a — [3z [16]. Note that the scaling 
function W(x) depends on the boundary condition. For the periodic boundary case, we had W(x) — f(x) where 
f(x) is the Airy distribution function in Eq. (2). For the free boundary condition, at least for the EW interface, 
W(x) — F(x) with F(x) being the F-Airy distribution function in Eq. (75). Note that it is not easy to guess the 
boundary conditions at the edges of a block of size £. Therefore, it is difficult to know the precise form of the scaling 
function W(x). However, it turns out that in the scaling regime (see later), only the large x behavior of W(x) matters. 
Now, both for the periodic (EW and KPZ) and the free boundary conditions (only EW), we have seen in Eqs. (57) 
and (77) that the scaling function W(x) has a Gaussian tail for large x, 



W(x) ~ exp [— bx 2 



(81) 



where the constant 6 = 6 for the periodic case and b = 3/2 for the free case. This suggests that the large x tail of 
the scaling function is generically Gaussian, though the constant b is nonuniversal and depends on the details of the 
boundary conditions. 

Let us first define the cumulative distribution of the global MRH in a system of size L and at time t, 



F(h m ,L,t) =Prob [m&x{h m (l),h m {2),...,h m (N)} < h m , L, t] . 
Using the fact that the block MRH /i m (i)'s have a factorised joint distribution, it then follows that 
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Substituting the scaling form, p(x,^) — £ a W(x^ a ) in Eq. (83) we get 



F(h m ,L,t) 





N 




W(y)dy 


i=a exp 








Jh 



W(y)dy 



(84) 



where the last equality holds near the tail h m — > oo. Using the Gaussian behavior of W(x) for large x in Eq. (81), it 
is easy to analyze the integral inside the exponential in Eq. (84) and one arrives at the limiting distribution, 



F(h m , L,t)^v[c^ (h m - D L )^j , with V(y) = exp [- e -»] , 



(85) 



where = £"-\/log(-L/£;), c is an unimportant nonuniversal constant and the scaling function V(y) is the celebrated 
Gumbel function. The probability distribution of the global MRH is then obtained by taking the derivative in Eq. 
(85), P(h m , L, t) = dF(h m , L, t)/dh m . This distribution is clearly peaked around its average value, 



(Mi)> ~ £Vlog(£/0 ~ ^ [log(£t~ 1/z ) 
The width of the peak around this average can also be read off Eq. (85) 

= V '<>&> (h m ) 2 ~ ^ ~ ^[log^r 1 /*) 
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(86) 



(87) 



For the (1 + l)-dimcnsional EW equation, z = 2, (3 = 1/4 and for the KPZ equation, z = 3/2, (3 = 1/3. The result 
for the average in Eq. (86) was also obtained in Rcf. [21], though the power 1/2 of the logarithm in Eq. (86) was 
found only numerically. In our case, this power 1/2 is a direct consequence of the large x Gaussian tail of the scaling 
function W(x) in Eq. (81), which were derived analytically for the EW and the KPZ equation in Section-IV and V. 
Thus, the main conclusion of this section is that in the growing regime, the appropriately scaled global MRH has the 
universal Gumbel distribution that does not depend on the details of the interface. 



VII. CONCLUSION 



To summarize, in this paper we have studied the distribution of the global maximum h m of relative heights (,i.e., 
the height measured with respect to the spatially averaged height) in (1 + l)-dimensional fluctuating interfaces, both 
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of the Edwards- Wilkinson and the Kardar-Parisi-Zhang variety. The distribution P(h m , L,t), at time t and in a 
system of size L, has two types of scaling behaviors depending on whether one is in the growing regime (where 
1 << t « L z ) or in the stationary regime (where t » L z ), z being the dynamical exponent of the interface. In 
the growing regime, the distribution (appropriately scaled) is described by the universal Gumbel law of the extreme 
value statistics. On the other hand, in the stationary regime, the distribution becomes time independent and has the 
scaling form, P(h m ,L) = L~ 1 / 2 W (/i m £~ 1//2 ) . The scaling function W(x) depends on the boundary condition. For 
the periodic boundary condition, we have analytically shown that W(x) — /(x) where /(x) is the Airy distribution 
function in Eq. (2) that describes the distribution of the area under a Brownian excursion over a unit interval. In 
fact, in this case, the EW and the KPZ interface both share the same scaling function f(x). In the case of the free 
boundary condition, we have shown that for the EW interface, W(x) = F(x) where F(x) is the -F-Airy distribution 
function defined in Eq. (75). These analytical results have also been verified by numerical integration of the EW and 
the KPZ equations. These results are summarized in Fig. 9. Since the relative heights in the stationary state are 
strongly correlated, our results provide a rather rare exactly solvable case for the distribution of extremum of a set of 
correlated random variables. 
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FIG. 9. As time t increases, the scaled MRH distribution crosses over from the universal Gumbel form in the growing regime 
to either the Airy distribution function (for the periodic boundary condition (PBC)) or to the F-Airy distribution function (for 
the free boundary condition (FBC)) in the stationary regime. 

Our work shows that there is yet another nontrivial scaling function in the (1 + l)-dimensional KPZ equation, 
namely the Airy distribution function that appears in the stationary MRH distribution in (1 + l)-dimensional KPZ 
equation with periodic boundary condition. Note that several other nontrivial scaling functions have appeared before 
in (1 + l)-dimensional KPZ equation. Let us list a few. 

(i) In the steady state of the one dimensional EW or the KPZ equation, the width of the surface height defined as 

2 

(88) 

was shown to have a scaling distribution, P(w 2 ) ~ (w 2 ) 1 ®(w 2 /(w 2 )) where the scaling function <&(x) is nontrivial 
[30]. 

(ii) For the (1 + l)-dimcnsional KPZ equation, the probability distribution of the spatially integrated height Y t = 
J H(x, t)dx, i.e., the area under the surface profile, was shown to have the following behavior at late times t [44] 

P(Y t ) ~ exp [tZ L (Y t /t)} , (89) 

where Zl(x) is a large deviation function which has a Gaussian peak near x = 0, but has highly non-Gaussian and 
asymmetric tails for large \x\. 

(hi) More recently, it has been shown [45] that the single site height distribution P(H,t) in the growing regime 
t « L 3 / 2 of the (1 + l)-dimensional KPZ equation scales as 

P{ H,t)~t-V*T w (^P-y (90) 

where the scaling function Tw (x) is identical (upto a scale factor) to the celebrated Tracy- Widom distribution for the 
largest eigenvalue of a random matrix drawn from the Gaussian unitary ensemble [46]. Note that in the stationary 
regime t » L 3 / 2 , this single site height distribution becomes a simple Gaussian. However, as demonstrated in this 
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paper, the global maximum of all these heights in a finite system of size L has a nontrivial distribution in the stationary 
regime, characterized by the Airy distribution function for the periodic boundary condition case. Determining the 
corresponding scaling function in the KPZ equation in a finite system with free boundary conditions remains an open 
problem. 

In this paper we have analytically studied the distribution of the extreme height fluctuations in the stationary state 
of a (1 + l)-dimensional interface. These height fluctuations are strongly correlated in the stationary state. Thus 
our work provides a rare example where one can calculate exactly the extreme value distribution of a set of strongly 
correlated random variables. In the context of interfaces, it is worth pointing out that another type of extreme 
statistics has been studied recently [47] . In Ref. [47] the authors study the distribution of the maximum of the Fourier 
modes of the height fluctuations. However, in this case the Fourier modes are independent of each other and the 
analytical computation is simpler, though physically relevant. 

Experimental Consequences: Our work opens up an interesting experimental challenge to measure the Airy (or the 
F-Airy) distribution function in a physical system. Many experimental systems are known to be well described by 
the (1 + l)-dimensional EW equation (28). Examples include, amongst others, the high-temperature step fluctuations 
in crystal surfaces [48,49] and the displacements of nonmagnetic particles in dipolar chains at low magnetic field 
[50]. In fact, in the former system, there have been beautiful recent measurements of a variety of interface properties 
such as various first-passage properties [48], which were studied before only theoretically [51]. Hence, it would be 
interesting to see if the MRH distribution can also be measured in these experimental systems. This would then be 
the first experimental measurement of the Airy distribution function which has so far appeared only in mathematical 
problems. 

Several new directions and interesting open questions emerge from this work. One of the challenges is to calculate 
the MRH distribution for interfaces in higher dimensions. For example, let us consider the (d + l)-dimensional EW 
equation, where the Laplacian term in the EW equation is (i-dimensional. This equation for d > 2 describes several 
interesting systems, e.g., a fluctuating random Gaussian manifold and also the time evolution of the magnetization 
of a spin system within mean field theory (for d > 4) . It is simpler to continue to think in terms of interface heights 
growing over a d-dimensional substrate. The relative heights will again reach a stationary joint distribution function 
in the long time limit, 



The calculation of the distribution of the maximal relative height h m , starting from the above joint distribution, is a 
theoretical challenge. For d < 2, the surface will be rough in the stationary state and one expects that the MRH will 
scale as, h m <~ L a where a = 1 — d/2 is the roughness exponent and L is the linear size of the substrate. Hence, the 
MRH distribution is expected to have the scaling form, P(h m ,L) <~ L~ a f d (h m L~ a ) where the scaling function fd(x) 
will depend explicitly on the dimension d. On the other hand, for d > 2, the surface is smooth and h m ~ O(l) even in 
the thermodynamic limit L — > oo. Hence, P(h m , L) will approach a limiting distribution P{h m ) in the thermodynamic 
limit. It would be interesting to compute this limiting distribution for d > 2. This is nontrivial because the relative 
heights are still correlated (the correlation function falls off only as a power law l/r d ~ 2 ), even though the surface is 
smooth. The only simplification happens in the limit of d — > oo, where the correlation function dampens exponentially 
and one would expect, using the theory of extreme value statistics of independent random variables as in Sec. VI, a 
Gumbcl distribution for P(h m ). In fact, this has recently been verified for the EW interface defined on a small world 
network [52]. However, for finite d > 1, the problem is wide open. The path integral techniques used in this paper are 
suitable only in one dimension. Entirely new techniques are required to deal with the problem in higher dimensions. 
Developments of such techniques are more than welcome. 
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APPENDIX A: DERIVATION OF TWO IDENTITIES 



In this appendix, we prove two identities 
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where a k 's are the magnitude of the zeros of the Airy function Ai(z) on the negative real axis, i.e., Ai(— a k ) = and 
B(a k ) is given by, 

f°° Ai(z)dz 

B ^ = Ah-a k) ■ ^ 

where Ai'(z) = dAi{z)/dz. The first identity was used in Eq. (24). The second identity appeared in Eq. (67) where 
C(a k ) = B\a k ). 

To prove the first identity in Eq. (Al), we consider the following contour integral in the complex z plane, 

r dz i STMz'W 

1 Jc 2m {z - u) Ai{z) ' 1 ' 

where the contour C is shown in Fig. 10. Note that the integrand in Eq. (A4) has simple poles at z = u and z = —a k 
for k = 1, 2, The latter fact follows since — a k is a simple zero of Ai(z). 



FIG. 10. The contour for the integral in Eq 
function Ai(z). 




The dots on the negative real axis denote the zeros — afe's of the Airy 



Using the asymptotic expansion of Ai(z) <~ z -1 / 4 exp[— 2z 2 / 3 /3] for large z [4], it is easy to see that the function 
Q(z) = Ai(z')dz' /Ai(z) <~ z -1 / 2 for large z. Hence, once we extend the radius of the contour C to infinity, the 
contribution to the integral I\ in Eq. (A4) from the circular part of the contour tends to zero. The only remaining 
contribution comes from the arms of the contour around the negative real axis, which is equal to the sum of the 
residues around — a k s each of which is a simple pole of the function Q(z). Evaluating the residues around these poles 
in Eq. (A4) we get 



B(a k ) 
(a k + u) ' 



(A5) 



On the other hand, since the integrand in Eq. (A4) is analytic in the interior of the domain between the contour 
C and Co in Fig. 10, it follows that J c = — J Co - But the integral around Co is just the residue at the pole z = u. 
Evaluating this residue, we then get 



= Su Mz)dz 
1 Ai(u) 

Comparing Eqs. (A5) and (A6), we get the identity in Eq. (Al). 

The identity in Eq. (A2) can be proved exactly in the same manner. We now consider the contour integral, 



(A6) 



r dz i 

J c 2ni (z - u) 



f?Ai{z') dz' 
Ai{z) 



-i 2 



(A7) 



around the same contour C in Fig. 10. Note that the integrand in Eq. (A7) now has doubles poles at z = —a k for 
all k = 1,2,... (apart from the single pole at z — u). The rest of the calculation is similar as in the evaluation of Ii 
and we do not repeat them here. Following these steps one easily arrives at the identity in Eq. (A2). 
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APPENDIX B: THE INVERSION OF A LAPLACE TRANSFORM 



In this appendix, we wish to express the function F(x) explicitly in real space x by inverting its Laplace transform, 

(Bl) 



f ' — ' 

/ F(x)e- sx dx = 2- 1 / 3 s 2 ^ VCK)e- 
Jo fc=i 



2/3,-1/3 



where the constant, C(a k ) = [B(afe)] with B{a k ) given in Eq. (A3). 

To proceed, we use an identity originally derived by Takacs in the context of Brownian excursions [3], 

2 1 ' 3 J J (2 1 / 3 x- 2 / 3 ) x- 5 / 3 e- sx dx = e- sV \ (B2) 

where the function J(x) can be expressed in terms of the known hypergeometric function U(a, b, z), 

2 2 ' 3 x 



J(x) = -f^tf (1/6,4/3, 2x 3 /27) e" 2 ^. 
Thus the identity in Eq. (B2) gives us the inverse Laplace transform, 
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2 1 / 3 a;' 5 / 3 j(2 1 / 3 a; - 2 / 3 ) 
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where J(x) is given by Eq. (B3). 

Our first step is to rescale s in Eq. (B4) by a constant factor, i.e., s — ► 7 3 / 2 s. This is equivalent to the rcscaling, 
x — > 7 _3 / 2 a;. Then, Eq. (B4) transforms into 
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Differentiating Eq. (B5) with respect to 7 gives 
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where J'(x) = dJ/dx. We take the derivative of the function J (a;) in Eq. (B3) and use the following identity [4] 
satisfied by the function U (a, b, x) 



dU(a,b,x) , , . TT . , . TT . j . 
x = (a — b + x) U (a, b,x) — U {a — 1, b, x). 



(B7) 



Rearranging and simplifying, we get 
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[U (1/2, 4/3, 4 7 3 /27a; 2 ) + 2U (-5/6, 4/3, 4 7 3 /27x 2 )] . (B8) 



Note that this formula in Eq. (B8) is exactly what we need to invert the Laplace transform in Eq. (Bl). Substituting 
7 = 2~ 1 / 3 a/ c in Eq. (B8) and using this formula in Eq. (Bl) we get 



2 -i/3 



F(x) = — V a k C(a k ) [U (1/6, 4/3, b k /x 2 ) + 2U (-5/6, 4/3, b k /x 2 )] e 
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where b k = 2a|/27. 
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